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A CRITICAL  EXAMINATION  OF  A NUMERICAL  FRACTURE  DYNAMIC  CODE 

by 

L.  Hodulak,  A.S.  Kobayashi  and  A.F.  Emery 

ABSTRACT 

After  upgrading  the  energy  dissipation  algorithm,  numerical  experi- 
ments were  conducted  to  assess  the  reliability  of  the  explicit  dynamic 
finite  element  code,  HCRACK.  Two  dynamic  fracture  specimens,  i.e.,  the 
wedge-loaded  rectangular  DCB  (RDCB)  specimen  and  the  wedge-loaded  tapered 
DCB  (TDCB)  specimen,  which  were  studied  experimentally  by  Kalthoff  et  al, 
were  then  analyzed  with  this  updated  fracture  dynamic  code.  Using  the 
experimentally  determined  dynamic  fracture  toughness,  KjD,  versus  crack 
velocity,  a,  relation,  the  RDCB  specimen  was  analyzed  first  by  the  "propa- 
gation method"  where  good  agreements  between  calculated  and  measured  KjD 
versus  a^  relations  were  observed.  The  calculated  a versus  time,  t, 
relation  was  then  used  as  input  data  in  the  "generation  method"  where  the 
resultant  KjQ  were  virtually  Identical  to  those  obtained  in  the  propaga- 
tion method.  Error  analyses  of  the  generation  method  were  also  made  first 
by  using  the  experimentally  determined  a versus  t relation  and  secondly  by 
artificially  perturbing  this  relation. 

A TDCB  specimen  was  then  analyzed  with  both  the  propagation  and  gen- 
eration methods  by  using  the  KjD  versus  a relation  established  for  this 
specimen  and  the  measured  a versus  £ relation,  respectively.  The 
computed  obtained  by  both  methods  were  in  good  agreement  with  the 

experimental  results,  showing  that  either  approach  can  be  used  In  analyzing 
fracture. 
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INTRODUCTION 

For  the  past  three  years,  two  of  the  authors  have  used  a two-dimension- 
al elasto-dynamlc  finite  element  code,  which  was  based  on  HONDO  [1],  to 
compute  the  dynamic  stress  Intensity  factor  for  a crack  propagating  with 
a prescribed  velocity  [2-5]  by  applying  to  each  node  a nodal  force 
sufficient  to  release  the  node.  This  numerical  procedure  was  later  modi- 
fied to  include  a startup  procedure  for  computing  dynamic  stress  intensity 
factor,  dynamic  energy  release  rate,  fracture  energy,  kinetic  energy  and 
strain  energy  at  each  increment  of  crack  advance  [6,7].  Also  the  impulse 
stress  waves  generated  by  the  instantaneous  application  of  a nodal  force 
to  model  the  release  of  a crack-tip  node  was  reduced  by  varying  the  force 
over  the  time  necessary  for  the  crack  tip  to  advance  one  nodal  distance. 
Physically,  this  procedure  models  a more  gradual  transit  of  the  crack-tip 
between  two  adjacent  finite  element  nodes  and  is  similar  to  that  developed 
by  Keegstra  [8-10]  with  the  exception  that  our  restraining  nodal  force  is 
completely  eliminated  when  the  crack-tip  reaches  its  adjacent  node.  Other 
nodal  force  release  mechanisms  Include  those  of  Malluck  and  King  [11]  and 
Rydholm,  Fredriksson  and  Nilsson  [12]  with  different  postulated  rates  of 
nodal  force  release.  The  dissipated  energy  during  a crack  extension  based 
on  any  of  the  above  three  nodal  force  release  mechanism  is  then  computed 
from  the  nodal  force  versus  nodal  displacement  relation  during  this  Incre- 
mental crack  extension.  In  general  this  nodal  force  versus  nodal  displace- 
ment relation  Is  non-linear  and  Is  goverened  by  the  dynamic  state  surround- 
ing the  propagating  crack  tip  thus  requiring  the  monitoring  of  nodal  force 
or  nodal  displacement  or  both  at  every  Incremental  time  in  the  dynamic 
finite  element  analysis.  Interestingly  enough,  recent  studies  showed  that 
the  differences  In  the  mechanism  of  nodal  force  release  [13,14]  caused 
little  changes  In  the  resultant  dynamic  stress  Intensity  factor.  It  Is 
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thus  of  no  surprise  that  good  to  excellent  agreements  were  claimed  by  all 
[6,11,12]  when  these  three  crack  tip  energy  dissipation  procedures  for 
computing  the  dynamic  stress  intensity  factor  was  checked  by  analyzing 
the  Broberg  problem  [15]. 

The  above  procedure  of  computing  dynamic  stress  intensity  factor  for 
a crack  whose  velocity  is  prescribed  to  be  equal  to  measured  one  in  a dy- 
namically fracturing  specimen  was  termed  "generation  calculation"  by 
Kanninen  [16,17]  who  also  expressed  reservations  on  the  accuracy  of  this 
numerical  approach.  The  "propagation  calculation"  in  contrast  to  the 
"generation  calculation"  is  based  on  an  assumed  dynamic  fracture  toughness, 
Kid,  versus  crack  velocity,  a,  relation  which  is  then  used  to  propagate 
a crack  [16-23].*  The  assumed  KjD  versus  a relation  is  considered  correct 
when  the  calculated  crack  propagation  history  coincides  with  the  experi- 
mental data,  and  the  «ID  at  crack  arrest,  if  any,  is  considered  to  be  the 
crack  arrest  toughness,  KJa,  sought  by  some  in  predicting  fracture  arrest 
of  a dynamically  propagating  crack. 

While  one  can  debate  the  merits  of  propagation  versus  generation  cal- 
culations, only  one  study  which  Involved  both  propagation  and  generation 
calculation  using  the  same  numerical  algorithm  [23]  has  been  published  to- 
date.  Since  the  limited  study  in  Reference  [23]  did  not  provide  a compre- 
hensive error  assessment  of  the  two  procedures,  this  paper  will  report 
on  our  comparative  studies  using  two  Araldite  B fracture  specimens  which 
were  analyzed  by  Kalthoff  et  al  by  the  method  of  caustics  [24,25]. 


* Note  that  Keegstra  In  References  [8,9,10]  used  the  propagation  calcu- 
lation to  compute  versus  a relations  in  fracturing  specimens. 


4 


DYNAMIC  FINITE  ELEMENT  ANALYSIS 

In  the  previous  studies  cited  above  [6,26],  the  dynamic  fracture  dy- 
namic code  HCRACK  was  shown  to  be  an  efficient  and  inexpensive  method  for 
simulating  dynamic  fracture  problems.  Numerical  experiments  proved  that 
reasonable  numerical  accuracy  can  be  obtained  by  using  coarse  meshes  of 
conventional  elements  (see  Figures  1 and  2)  and  a moderate  number  of  time 
steps,  e.g.,  about  150  steps  for  crack  propagation  and  subsequent  arrest 
in  a RDCB  specimen  shown  In  Figure  1.  Unlike  the  implicit  dynamic  finite 
element  codes  used  by  others,  however,  it  was  difficult  to  accurately  pre- 
scribe the  rate  of  nodal  force  release  since  the  input  nodal  force  would 
not  generally  be  in  equilibrium  with  the  dynamic  state  of  stress  in  the 
adjoining  finite  elements  in  this  explicit  finite  element  code.  As  a re- 
sult, an  in-depth  study  on  the  performance  of  our  fracture  dynamic  code 
was  conducted  for  different  crack  tip  force  release  rates,  different  cal- 
culation procedures  for  the  dynamic  stress  intensity  factor,  and  different 
finite  element  breakdown.  A brief  description  of  some  of  these  findings 
are  presented  in  the  following. 

As  mentioned  above,  the  algorithm  for  artificially  prescribing  an  in- 
put nodal  force  at  the  crack  tip  for  each  time  step  for  prescribed  decrease 
in  the  resultant  residual  nodal  force  in  the  dynamic  code  is  not  straight- 
forward and  often  a complete  release  of  the  nodal  force  cannot  be  achieved 
in  the  prescribed  time  period.  The  basis  of  the  numerical  method  is  to 
define  the  force,  F 1 which  must  be  applied  to  a node  at  time  step  n such 
that  the  time  variation  of  the  stress  follows  the  form  shown  on  figure  3. 

In  an  Implicit  code,  application  of  Fn  would  result  In  the  same  calculated 
force  at  the  end  of  the  time  step.  With  the  explicit  code,  however,  the 
calculated  force  at  the  end  of  the  Increment,  Fn°,  will  rarely  be  equal  to 
F *.  Accordingly,  the  force  at  the  next  time  step  Fnll  must  be  adjusted 


5 


not  only  to  compensate  for  the  error,  but  also  to  yield  the  desired  value 
at  the  end  of  the  time  step.  The  simplified  method  used  in  [26]  was  re- 
placed by  the  following  equation  and  typical  results  are  as  shown  in 
figure  3. 


n F 1 

i _ c 0 _ c prescribed  _ j 

n+1  n n+l  " . ^ 2n"J+^ 


0) 


At  the  beginning  of  the  crack  propagation  history  when  the  first  crack 
tip  node  is  relased  (see  Figure  3),  excellent  agreement  between  the  pre- 
scribed (linear  decrease  in  this  case)  and  actually  achieved  nodal  forces 
is  noted  while  for  the  sixth  node  which  was  released  later  (see  Figure  3), 
some  deviations  between  both  nodal  forces  are  noted.  Initial  static 
equilibrium  prior  to  crack  propagation  most  likely  contribute  the  excellent 
results  in  the  former  case. 

Also  noteworthy  is  the  recent  study  by  Malluck  and  King  [13]  who 
compared  energy  release  rates  for  the  two  distinctly  different  functions 
of  F/Fq  = [1-b/A]3^2  and  F/Fq  = [1-b/A]^2,  where  b is  the  distance  between 
hypothetical  crack  tip  location  and  the  released  crack  tip  node  and  A is 
the  inter-nodal  distance  and  F and  Fq  are  the  instantaneous  and  original 
crack  tip  nodal  forces,  respectively.  Their  results  showed  no  significant 
differences  in  the  calculated  dynamic  stress  intensity  factors  for  crack 
speeds  lower  than  25  percent  of  the  shear  wave  velocity,  i.e.  c < .25  C£. 
Our  use  of  a linearly  decreasing  nodal  force,  F/Fq  = [T-b/A],  with  constant 
crack  velocity  between  the  two  adjoining  finite  element  nodes  is  thus 
justified. 

The  dynamic  stress  Intensity  factor  was  computed  directly  by  the 
total  strain  energy  released  from  an  instantaneous  balance  of  the  total 
energy  of  the  entire  specimen  [7]  as 


where  En‘  En+1  are  the  total  strain  energies  for  crack  lengths  of  an, 
an+i’  respectively  when  the  crack  extended  from  node  n to  node  n+1.  The 
dynamic  stress  intensity  factor,  Kj,  was  then  computed  from  Gj  using 
Freund's  relation  [24].  Alternatively  the  value  of  Gj  was  computed  by 
energy  dissipated  at  the  released  node  as 

m 

Gj  = (u.  A Fi  + (u1  - 4F()/  (an+1  - an)  (3) 

where  m is  the  number  of  time  steps  between  nodes  n and  n+1,  u^  and  a F^ 
are  displacement  and  decrease  of  force  at  the  released  node  n,  respectively. 

Figure  4,  shows  the  dynamic  fracture  toughness,  Kj0,  associated  with 
crack  propagation  and  arrest  in  one  of  Kalthoff's  RDCB  specimens  [24]  com- 
puted by  both  equations  (2)  and  (3)  using  the  "propagation  method." 

Although  details  of  this  analysis  are  described  in  the  following  section, 
the  results  are  shown  in  this  section  as  an  indication  that  little  differ- 
ence can  be  noted  in  the  KID  obtained  by  the  two  algorithms. 

As  shown  in  Figure  4,  the  forced  linear  decrease  in  the  crack  tip 
nodal  force  improved  the  simulation  of  the  smoothly  propagating  crack  and 
eliminated  the  spurious  oscillations  in  dynamic  stress  intensity  factor 
observed  previously  [2-5].  It  Is  uncertain,  however,  to  what  extent  this 
smoothing  procedure  may  hide  the  true  oscillations  of  the  dynamic  stress 
intensity  factor  eventually  Induced  by  the  reflected  stress  waves  which 
emanated  from  the  running  crack. 

SPECIMENS  AND  MATERIAL  DATA 

The  two  specimens  analyzed  by  the  dynamic  finite  element  code  are 
the  welge-loaded,  RDCB  and  TDCB  specimens  which  were  investigated  experi- 
mentally by  Kalthoff  et  al  [24,25].  Specimen  geometries  of  these 
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Araldite  B specimens  and  their  finite  element  idealizations  can  be  seen 
in  Figures  1 and  2. 

Although  the  rigid  loading  wedge  between  the  two  loading  pins  will 
prevent  any  inward  displacement  of  the  loading  pins,  these  pins  are  free 
to  leave  the  wedge  and  travel  outwards.  The  resultant  dynamic  stress  in- 
tensity factors  in  the  presence  of  separating  pins  could  vary  significant- 
ly during  crack  propagation  [26].  The  smaller  mass  density  and  the  two 
orders  of  magnitude  larger  compliance  of  the  Araldite  B specimens  in 
comparison  to  the  steel  specimen  studied  in  Reference  [26]  should  have 
reduced  the  additional  input  energy  due  to  any  possible  separation  of  the 
loading  pins  and  thus  constant  loading  pin  displacement  were  prescribed 
at  the  pin  holes. 

Material  constants  of  Araldite  B used  for  this  dynamic  finite  element 
analysis  after  [24]  are  modulus  of  elasticity  E = 3.38  GP  . Poisson  ratio 

a 

of  v = 0.33  and  mass  density,  p = 1047  kg/m3.  The  experimentally  determined 
dynamic  fracture  toughness  K^,  versus  crack  velocity,  a,  relations  used 
in  the  propagation  calculations  of  RDCB  and  TDCB  specimens  are  both  plotted 
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Kalthoff's  measured  oscillating  Kjp  values  [25]  after  crack  arrest,  was 
used  in  the  analysis  of  the  TDCB  specimen. 

RESULTS 
RDCB  Specimen 

The  first  numerical  analysis  involved  a propagation  calculation  for 
the  RDCB  specimen  of  Figure  1 using  the  KjQ  versus  a relation  of  Figure  5 
and  a Kjq  = 2.32  MNm  . The  resulting  dynamic  fracture  toughness  and 
crack  tip  motion  of  this  propagation  calculation  are  shown  in  Figures  6 
and  7,  respectively.  The  "propagation"  crack  tip  motion  from  Figure  7 was 
then  used  as  input  data  for  the  "generation"  calculation.  This  result  is 
not  plotted  in  Figure  6 since  the  Kjp  versus  k relations  obtained  by  both 
the  propagati™  and  generation  calculation  were  indistinguishable. 

As  an  additional  numerical  experimentation,  however,  the  measured 
crack  length,  a,  versus  time,  t,  relation  of  Reference  [24]  was  used  as 
input  to  the  "generation"  calculation  and  the  resultant  KjD  are  also  shown 
in  Figures  6 and  7.  Despite  the  lack  of  complete  agreement  between  the 
two  Kjp  curves  obtained  by  propagation  and  generation  calculation,  the 
shapes  of  these  two  curves  are  very  close.  Although  both  Kjp  curves  agree 
well  with  experimental  data  during  the  first  half  of  dynamic  crack  propa- 
gation as  shown  in  Figure  7,  a distinct  difference  is  noted  by  a second 
local  maximum,  which  occurs  In  both  propagation  and  generation  calculations 
prior  to  crack  arrest,  but  which  does  not  occur  in  the  experimental  results. 
The  similarity  between  the  propagation  and  generation  KjD  curves  Is  more 
apparent  In  Figure  8 where  the  second  maxima  In  the  two  calculations  occur 
at  the  same  time.  The  higher  KjD  values  In  the  generation  calculation  at 
lower  measured  crack  velocities  during  much  of  the  crack  propagation  will  re- 
sult in  a general  shift  of  two  KjQ  versus  a relation  in  Figure  5. 
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Figure  6 also  shows  that  the  computed  crack  jump  distance  is  4% 
shorter  of  the  measured  one  in  the  propagation  calculation  but  by  defini- 
tion is  equal  to  measured  distance  in  the  generation  calculation.  Although 
the  propagation  calculation  is  terminated  when  the  computed  dynamic  stress 
intensity  factor  falls  below  the  minimum  Kjq  value  in  Figure  5,  the  gener- 
ation calculation  is  continued  up  to  the  prescribed  crack  tip  length  and 
crack  arrest  time.  Significantly  lower  dynamic  stress  intensity  at  the 
instant  of  crack  arrest  is  noted. 

The  sensitivity  of  the  dynamic  stress  intensity  factor,  which  is  cal- 


culated by  the  generation  method,  to  the  instantaneous  crack  velocity  is 
further  demonstrated  in  Figure  8.  In  order  to  assess  the  sensitivity  of 
Kid  obtained  by  the  generation  method  to  the  input  data,  a numerical  experi- 
ment was  conducted  by  artificially  perturbing  the  smooth  experimental 
curve  of  the  crack  tip  motion  in  Figure  8.  The  result  was  a severely 
perturbed  KjQ  also  shown  in  Figure  9,  where  discrete  increases  and  de- 
creases in  crack  velocities  are  followed  by  local  minima  and  maxima  of 
Kid  respectively. 

TDCB  Specimen 

Figure  9 shows  the  KjD  as  a function  of  a computed  by  the  propagation 
method,  using  the  Kjg  versus  a relation  of  Figure  5 and  by  the  generation 
method  using  experimentally  determined  a versus  t relations  for  the  TDCB 
specimen  together  with  experimental  data  from  Reference  [25].  A second 
maximum,  which  resembles  that  found  previously  in  the  RDCB  specimen,  in 
Kid  can  be  observed.  The  computed  crack  jump  distance  obtained  by  the 
propagation  method  is  shorter  than  the  experimental  one  by  12%.  In  the 
propagation  calculation  the  computed  K^q  Increased  again  to  a value 
approaching  experimental  Kjp  after  the  Initial  crack  arrest.  Subsequently 
computed  KID  oscillated  around  the  few  experimental  points. 


Figure  10  shows  the  Kjg  versus  t relations  obtained  by  both  propagation 
and  generation  calculations.  Although  the  two  calculated  Kjq  are  in 
excellent  agreement  with  each  other  except  for  the  initial  phase  of  crack 
propagation  in  this  TDCB  specimen,  the  calculated  KID  are  lower  than  the 
measured  KjQ  just  prior  to  and  after  crack  arrest.  Previous  experience 
with  steel  TDCB  specimens  [4,5,26]  indicate  that  this  small  underestimate 
could  be  attributed  to  the  possible  separation  of  the  loading  pins  from 
the  loading  wedge  during  crack  propagation. 

CONCLUSIONS 

The  results  of  the  present  and  of  the  previous  studies  using  HONDO  II 
show  that  the  dynamic  stress  intensity  factor  for  a crack  propagating  in 
a finite  two-dimensional  body  can  be  computed  relatively  inexpensively  with 
an  accuracy  sufficient  for  many  practical  purposes.  Very  close  agreements 
between  the  KID  obtained  by  the  generation  and  by  the  propagation  calcula- 
tions should  dispel  the  reservations  [16,17]  about  this  dynamic  fracture 
algorithm. 

When  used  in  conjunction  with  measured  crack  position  versus  time 
data,  the  generation  method  with  proper  care  can  be  used  to  accurately 
calculate  the  dynamic  stress  intensity  factor  during  the  fast  crack  propa- 
gation and  crack  arrest. 

On  the  other  hand  the  uncertainty  in  the  KjQ  versus  a relations, 
particularly  in  the  region  of  very  low  velocities  together  with  limitation 
in  the  finite  element  modeling  of  dynamic  crack  propagation  offers  little 
chance  for  simulating  the  crack  propagation  and  crack  arrest  event  by  the 
propagation  method  when  the  dynamic  stress  intensity  factor  oscillates  in 
a narrow  range  about  the  crack  arrest  stress  intensity  factor  as  shown  by  some 
experimental  results  with  the  single  edged  notch  specimens  reported  in  [25]. 
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DISCUSSION 

It  has  been  a common  practice  by  all,  including  the  authors,  to 
verify  their  fracture  dynamic  code  by  analyzing  the  Broberg  problem  [15] 
for  which  the  dynamic  solution  is  available.  Good  agreements  in  these 
studies  cannot  be  construed  as  verification  of  numerical  solutions  gener- 
ated for  cracks  propagating  in  finite  specimens  composed  of  real  materials. 
The  discrepancies  between  the  computed  and  the  experimentally  determined 
Kjg-values  shown  in  Figures  6 and  9 could  have  arisen  from  the  viscous 
damping  in  Araldite  B which  was  not  modeled  in  the  elasto-dynamic  analyses 
described  in  this  paper.  A study  of  the  time-dependent  energy  balance 
during  crack  propagation  and  arrest  suggests  that  the  consistently  appear- 
ing second  maxima  in  the  calculated  K^-curves  are  real  phenomena  based  on 
elastic  analyses.  It  is  interesting  to  note  that  the  limited  experimental 
KI0  versus  a relation  obtained  for  RDCB  specimens  machined  from  high 
strength  steel  [25]  is  in  qualitative  agreement  with  our  elastic  analysis 
of  the  RDCB  specimen. 
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FIGURE  I.  WEDGE- LOADED  RDCB  SPECIMEN. 
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FIGURE  2.  WEDGE-LOADED  TDCB  SPECIMEN. 
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FIGURE  6. DYNAMIC  FRACTURE  TOUGHNESS  IN  WEDGE-LOADED 
RDCB  SPECIMEN. 
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